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Abstract. We introduce a diffuse interface model for the phenomenon of electrowetting on di- 
electric and present an analysis of the arising system of equations. Moreover, we study discretiza- 
tion techniques for the problem. The model takes into account different material parameters on 
each phase and incorporates the most important physical processes, such as incomprcssibility, 
electrostatics and dynamic contact lines; necessary to properly reflect the relevant phenomena. 
The arising nonlinear system couples the variable density incompressible Navier-Stokes equations 
for velocity and pressure with a Cahn-Hilliard type equation for the phase variable and chemical 
potential, a convection diffusion equation for the electric charges and a Poisson equation for the 
electric potential. Numerical experiments are presented, which illustrate the wide range of effects 
the model is able to capture, such as splitting and coalescence of droplets. 



1. Introduction 

The term electrowetting on dielectric refers to the local modification of the surface tension 
between two immiscible fluids via electric actuation. This allows for change of shape and wetting 
behavior of a the two-fluid system and, thus, for its manipulation and control. 

The existence of such a phenomenon was originally discovered by Lippmann |39j . more than 
a century ago (see also [5J H3J M 157]). However, only recently has electrowetting found a wide 
spectrum of applications, specially in the realm of micro- fluidics |15l 1161 [25] . One can mention, for 
example, reprogrammable lab-on-chip systems [371 [55], auto- focus cell phone lenses [TU], colored oil 
pixels and video speed smart paper j33J [501 [5T] . In [35] , the reverse electrowetting process has been 
proposed as an approach to energy harvesting. 

From the examples presented above, it becomes clear that it is very important for applications to 
have a better understanding of this phenomenon and it is necessary to obtain reliable computational 
tools for the simulation and control of these effects. The computational models must be complete 
enough, so that they can reproduce the most important physical effects, yet sufficiently simple that 
it is possible to extract from them meaningful information in a reasonable amount of computing 
time. Several works have been concerned with the modeling of electrowetting. The approaches 
include experimental relations and scaling laws [34, 62 , empirical models |41j . studies concerning 
the dependence of the contact angle ([2H[55]) or the shape of the droplet ([1H[T7]) on the applied 
voltage, lattice Boltzmann methods [U [3J and others. Of relevance to our present discussion are 
the works [5H [S3] and [201 I23J- To the best of our knowledge, [SH [S3] are the first papers where 
the contact line pinning was included in an electrowetting model. On the other hand the models 
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of [201 I23| are the only ones that are intrinsically three dimensional and do not assume any special 
geometric configuration. They have the limitation, however, that they assume the density of the 
two fluids to be constant and they apply a no-slip boundary condition to the fluid-solid interface, 
thus limiting the movement of the droplet. 

The purpose of this work is to propose and analyze an electrowetting model that is intrinsically 
three-dimensional; it takes into account that all material parameters are different in each one of 
the fluids; and it is derived (as long as this is possible) from physical principles. To do so, we 
extend the diffuse interface model of [20]. The main additions are the fact that we allow the 
fluids to have different densities - thus leading to a variable density Cahn Hilliard Navier Stokes 
system - and that we treat the contact line movement in a thermodynamically consistent way, 
namely using the so-called generalized Navier boundary condition (see [JHIHH])- In addition, we 
propose a (phenomenological) approach to contact line pinning and study stability and convergence 
of discretization techniques. In this respect, our work also differs from [20, 23], since our approach 
deals with a practical fully discrete scheme, for which we derive a priori estimates and convergence 
results. 

Through private communication we have become aware of the following recent contributions: 
discretization schemes for the model proposed in [2U] are studied in [35] ; the models of [20] H3] have 
been extended, using the techniques of [I], in [I!|1[2H] where discretization issues are also discussed. 

This work is organized as follows. In §1.1 1 we introduce the notation and some preliminary 
assumptions necessary for our discussion. Section [2] describes the model that we shall be concerned 
with and its physical derivation. A formal energy estimate and a formal weak formulation of our 
problem is shown in section [3] The energy estimate shown in this section serves as a basis for 
the precise definition of our notion of solution and the proof of its existence. The details of this 
are accounted for in section [4] In section [5] we discuss discretization techniques for our problem 
and present some numerical experiments aimed at showing the capabilities of our model: droplet 
splitting and coalescence as well as contact line movement. Finally, in section [6] we briefly discuss 
convergence of the discrete solutions to solutions of a semi-discrete problem. 




Figure 1.1. The basic configuration of an electrowetting on dielectric device [HI 
[T6] . The solid black region depicts the dielectric plates and the white region denotes 
a droplet of one fluid (say water) , which is surrounded by another (air) . We denote 
by fl the fluid domain, by T its boundary, by fi* the region occupied by the fluids 
and the plates and by d*Q* :— dil* \ T. 

1.1. Notation and Preliminaries. Figure [O] shows the basic configuration for the electrowetting 
on dielectric problem. We use the symbol to denote the domain occupied by the fluid and f2* for 
the fluid and dielectric plates, thus, fi C f2*. In this manner, we assume that f2 and fi* are convex, 
bounded connected domains in R d , for d = 2 or 3, with C 0,1 boundaries. The boundary of ft is 
denoted by T and d*£l* = dfl* \ T, n stands for the outer unit normal to T. We denote by [0,T] 
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with < T < oo the time interval of interest. For any vector valued function w : — > M. d that is 
smooth enough so as to have a trace on T, we define the tangential component of w as 

(1.1) w T \ r := w| r - (w| r -n)n, 

and, for any scalar function /, d T f :— (V/) x . 

We will use standard notation for spaces of Lebesgue integrable functions L p (il), 1 < p < oo and 
Sobolev spaces W™(S1) 1 < p < oo, to 6 No, [2] • Vector valued functions and spaces of vector valued 
functions will be denoted by boldface characters. For S C M d , by (-,-)s we denote, indistinctly, 
the L 2 (S)- or L 2 (S')-inner product. If no subscript is given, we assume that the domain is Q. If 
S C then the inner product is denoted by [•, -]s and if no subindex is given, the domain must 

be understood to be T. We define the following spaces: 

(1.2) := {v G H 1 ^*) : v\ d * Q * = 0} , 
normed by 

\MhI ■= HVuHi^n*), 

and 

(1.3) V := {v G H x (0) : vn| r = 0} , 
which we endow with the norm 

Mv ■■= l|V V ||2 3 + |K||^ 3(r) . 

Clearly, for these norms, they are Hilbert spaces. 

To take into account the fact that our problem will be time dependent we introduce the following 
notation. Let E be a normed space with norm || • \\e- The space of functions <p : [0, T] —¥ E such 
that the map (0,T) 9 t i-> \\<p{t)\\ E G R is LP-integrable is denoted by LP(0,T,E) or L P (E). To 
discuss the time discretization of our problem, we introduce a time-step At > (for simplicity 
assumed constant) and let t n — nAt for < n < N := \T/At\. for any time-dependent function, 
tp, we denote <p n :— (p(t n ) and the sequence of values {<p n }n=o * s denoted by <pAt- For any sequence 
fAt we define the time-increment operator by 

(1.4) <)ip n := ip n — <£ n-1 , 
and the time average operator (•) by 

(1.5) W-=\W n + ^ n - 1 )- 
On sequences <pAt C E we define the norms 

N N 

\\VAt\\% iE ) ~ A *5Z H^lll' IIPAt||*»(E) := max . {|^ n || E }, ||¥>At||^i/a (E) :=^||0^ 



0<n<Ar 

n— n— 1 



n||2 



which are, respectively, discrete analogues of the L 2 (E), L°°(E) and H 1 ^ 2 (E) norms. When dealing 
with energy estimates of time discrete problems, we will make, without explicit mention, repeated 
use of the following elementary identity 

(1.6) 2a(a-b) = a 2 - b 2 + {a - b) 2 . 
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2. Model Derivation 

In this section we briefly describe the derivation of our model. The procedure used to obtain 
it is quite similar to the arguments used in [20l |49l [T] and it fits into the general ideological 
framework of so-called phase-field models. In phase-field methods, sharp interfaces are replaced by 
thin transitional layers where the interfacial forces are now smoothly distributed and, thus, there 
is no need to explicitly track interfaces. 

2.1. Diffuse Interface Model. To develop a phase-field model, we begin by introducing a so- 
called phase field variable (f> and an interface thickness 6. The phase field variable acts as a marker 
that will be almost constant (in our case ±1) in the bulk regions, and will smoothly transition 
between these values in an interfacial region of thickness S. Having introduced the phase field, all 
the material properties that depend on the phase are slave variables and defined as 

(2.1) m arctan ' + * 2 

where the \&j are the values on each one of the phases. 



Remark 2.2 (Material properties). Relation (2.1| is not the only possible definition of the phase 



dependent quantities. For instance, 60J proposes to use a linear average between the bulk values. 
This approach has the advantage that the derivative of a phase-dependent field with respect to 
the phase (expressions that contain such quantities appear repeatedly) is constant, which greatly 
simplifies the calculations. However, this definition cannot be guaranteed to stay in the physical 
range of values which might lead to, say, a vanishing density or viscosity. On the other hand, [ID] 
proposes to use a harmonic average which guarantees that positive quantities stay bounded away 



from zero. In this work, we will assume that, with the exception of the permittivity e, (2.1) is the 
way the slave variables are defined, which has the advantage that guarantees that the field stays 
within the physical bounds. Any other definition with this property is equally suitable for our 
purposes. 

We model the droplet and surrounding medium as an incompressible Newtonian viscous two- 
phase fluid, so that its behavior is governed by the variable density incompressible Navier Stokes 
equations. The equation of conservation of momentum can be written in several forms. We chose 
the one proposed by Guermond and Quartapelle ( 30J, see also [551 150] ) because its nonlinear term 
possesses a skew symmetry property similar to the constant density Navier Stokes equations, 

(2.3a) cr(cru) t + (pu-Vu) + ~V(pu)u - V (??S(u)) + Vp = F, 

(2.3b) Vu = 0, 

where a = *J~p and p is the density of the fluid and depends on the phase field; u is the velocity of 
the fluid; p is the pressure; r\ is the viscosity of the fluid and depends on <j>; S(u) = |(Vu + Vu T ) is 
the symmetric part of the gradient and F are the external forces acting on the fluid. 

The phase field can be thought of as a scalar that is convected by the flow. Hence its motion is 
described by 

(2.4) & + V-fau) = -VJ*, 

for some flux field which will be found later. 

To model the interaction between the applied voltage and the fluid we introduce the charge 
density q. Another possibility, not explored here, is to introduce ion concentrations, thus leading 
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to a Nernst Planck Poisson-like system, see [23J HSU ESI H5] • The electric displacement field D is 
denned in CI*. The evolution of these two quantities is governed by Maxwell's equations, i.e., 

(2.5) V-D = g, D t + gu + J D =0, 

for some flux Jd ■ Notice that we assume that the magnitude of the velocity of the fluid is negligible 
in comparison with the speed of light, and that the frequency of voltage actuation is sufficiently 
small, so that magnetic effects can be ignored. Taking the time derivative of the first equation and 
substituting in the second we obtain 

(2.6) & + V(gu) = -VJ D . 

To close the system, we must prescribe boundary conditions, determine the force F exerted on 
the fluid, and find constitutive relations for the fluxes and Jd- We are assuming the solid walls 
are impermeable, therefore if n is the normal to L, un = on T and J[,-n = for any flux J[,. 
To find the rest of the boundary conditions, F and relations for the fluxes, we denote the surface 
tension between the two phases by 7 and define the Ginzburg-Landau double well potential by 

f(£ + i) 2 , £<-i, 

mo = j *(w 2 ) 2 , iei<i, 

U-1) 2 , £>i. 

Remark 2.7 (The Ginzburg Landau potential). The original definition, given by Cahn and Hilliard, 
of the potential is logarithmic. Sec, for instance, [27]. This way, the potential becomes infinite if 
the phase field variable is out of the range [—1,1], thus guaranteeing that the phase field variable 
<fi stays within that range. This is difficult to treat both in the analysis and numerics and hence 
practitioners have used the Ginzburg-Landau potential c(l — £ 2 ) 2 , for some c > 0. We go one step 
further and restrict the growth of the potential to quadratic away from the range of interest. With 
this restriction Caffarelli and Miiller, [13], have shown uniform L°°-bounds on the solutions of the 
Cahn Hilliard equations (which as we will see below the phase field must satisfy). This has also 
proved useful in the numerical discretization of the Cahn Hilliard and Cahn Hilliard Navier Stokes 
equations, see [5^1 [551 [55] . 

Finally, we introduce the interface energy density function, which describes the energy due to 
the fluid-solid interaction. Let 9 S be the contact angle that, at equilibrium, the interface between 
the two fluids makes with respect to the solid walls (see [1H1 HSJ [S3]) and define 

C0S6» S . / 7T<A 

e /s (0) = — sm^ T j. 

Then, up to a constant, the interfacial energy density equals 70/ s ((/>). 
Let us write the free energy of the system 

m *=-,J a (iw + iw W ) + 7 / r e,. ( « + \j a ^idi* + \j a ^u + \ I/. 

where e is the electric permittivity of the medium and A > is a regularization parameter. Com- 
puting the variation of the energy (£ with respect to 0, while keeping all the other arguments fixed, 
we obtain that 

(D^J) = f ^+ [ L4>, 
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where p is the so-called chemical potential which, in this situation, is given by 

(2.9) p = 7 Qw'W - SAc^j - ^|^|D| 2 + Ip'(<P)\u\ 2 . 
The quantity L is given by 

(2.10) L = j(e' f M + 5d n( f>), 
and can be regarded as a "chemical potential" on the boundary. 

Remark 2.11 (Chemical potential). From the definition of the chemical potential p we see that 
the product /xV0 includes the usual terms that define the surface tension, i.e., 



Additionally, it has the term 



- ( iw'(0)-(SA0) V-, 



2e(0) 2 ' 

which, in some sense, can be thought of as coming from the Maxwell stress tensor. 

With this notation, let us take the time derivative of the free energy: 
d£ 



/ iuj) t + / Lct) t + / E-D,+ / p(0)u-u f + A / qqt, 
Jn Jr Jn Jn Jn 



dt 

where E is the electric field, defined as E := e~ 1 D. Let us rewrite each one of the terms in this 
expression. Using (2.4) and the impermeability conditions, 

/ n4>t = - / (iV- (</>u + J^) = / Wp- ((j>u + 3$) . 
Jo Jo Jn 



Using (2.51 



/ E D t =- f E-{qu + J D ). 

Jn Jn 



For the boundary term, we introduce the material derivative at the boundary <f> = <f> t + UT-d T (j) and 
rewrite 

J L<j) t = J L(<p - \i T d T 4>). 
Notice that a(au) t — pu t + \ptU-, so that using (2.3), and integrating by parts, we obtain 

f p(4>)u-u t = /p.u-i / p{4>)M 2 + [ »j(S(u).n).u T - f 7 ? |S(u)| 2 . 
Jn Jn z Jn Jr Jn 

Finally, using ( |2.6[ ) and the impermeability condition (qu + Jrj)-n|r = 0, 

A / qq t = -A / gV- (gu + J D ) = A / Vq- (qu + J D ) . 

JO JO JO 

With the help of these calculations, we find that the time-derivative of the free energy can be 
rewritten as 



(2.12) 



<£=- / /iV0-u+ / J -Vm+ / L(4> - u T d T <t>) - / E-(gu + J D )+ / F-u 

Jo Jo Jr Jo Jo 

-~ / p'm t \u\ 2 + I *7(S(u)-n).u T - / 77|S(u)| 2 + ^ / u-V (q 2 ) + A / V<z-J D . 
* Jo Jr Jo z Jn Jn 
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From (2.12), we can identify the power of the system, i.e., the time derivative of the work 2U of 
internal forces, upon collecting all terms having a scalar product with the velocity u, 

227 = / Fu- [ /iV^-u- f ql]-u + ^V{q 2 )-u- \ [ p'^fau-u. 
Jn Jn Jn ^ * Jn 

We assume that the system is closed, i.e., there are no external forces. This implies that 211=0 
and we obtain an expression for the forces F acting on the fluid, 

F = pS7<t> + qE + i//($&u - V Qg 2 

Using the first law of thermodynamics 

d€ _ d2U d© 
"dT ~~ ~df ~ ~dT' 

where the absolute temperature is denoted by T and the entropy by 6, we can conclude that 

T&= [ ?/|S(u)| 2 - f E-J D + f J -V M + A / Vg-J D + / 77(S(u)-n)-u T + f L(<j) - u T d T <j>). 
Jn Jn Jn Jn Jr Jr 

To find an expression for the fluxes we introduce, in the spirit of Onsager [44, 49J, a dissipation 
function $. Since this must be a positive definite function on the fluxes, the simplest possible 
expression for a dissipation function is quadratic and diagonal in the fluxes, e.g. 

<S>=1 f ilJJ 2 + ? fi 2 + lf 4lJnl 2 + ^ f B\u ' 2 



2./o M rfl ' 2 ./ r r ' 2./o K> D| 2 



where all the proportionality constants, in principle, can depend on the phase <f>. Here, M is known 
as the mobility, K the conductivity and j3 the slip coefficient. Using Onsager's relation 

£j(<S(J) + $(J)Vj) = 0, vj, 



and (2.12), we find that 

(2.13) J = -MVp, J = -MV/i, J D = K (E - Wq) , /3u T = -7 ? S(u) nx + Ld T <f>, a<j) = -L, 
where S(u) nx := (S(u)-n) T . 



Remark 2.14 (Constitutive relations). Definitions (2.13) can also be obtained by simply saying 
that the constitutive relations of the fluxes depend linearly on the gradients, which is implicitly 
postulated in the form of the dissipation function $. 

Since, in practical settings, there is an externally applied voltage (which is going to act as the 
control mechanism) we introduce a potential V and then the electric field is given by E = — VU 
with V = Vq on <9*f2* , where Vq is the voltage applied. 

To summarize, we obtain the following system of equations for the phase variable <f> and the 
chemical potential p, 

{cj>t + u-V^ = V-(M(0)Vju), in O, 

p = 7 QvV(0) - 8A<fr) - \e'{4>)\^V\ 2 + ^'(^)|u| 2 , in O, 
a {fa + u T d T <f>) + 7 (Q'fsW + Sdnt) = 0, M(<f>)d n p = 0, on T, 
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and the velocity u and pressure p, 



(2.16) 



[ D(p(£)u) _ y (??WS(u)) + Vp = ^ -qV{V + Xq) + \p'{<t>)4>tu on 0, 



V-u = 0, 
un = 0, 

[P(<P)u T + r/(0)S(u) nT = 7 (0' f M + 5d n 4>) 8 T ^ 



in Q, 

on r, 
on r, 



where we have set 



D{P ^ )U) := ( r(0)(a(0)u) t + p^u-Wu + |V- W)u)u. 



In addition, we have the equation for the electric charges q, 



(2.17) 



f ft + V-(gu) = V [K(<t>)V {Xq + V)] , in 0, 
jiq»V(Ag + V)-n = 0, on T, 



and voltage V, 



(2.18) 



-V- (e*(0)VV) 

= 0, 



in Q*, 
on 9*0*, 

on dn* n r, 



where 



'e(0), n, 



with en being the value of the permittivity on the dielectric plates Q* \ tt, so Ed is constant there. 



Remark 2.19 (Generalized Navier Boundary condition). In (2.16), the boundary condition for the 
tangential velocity is known as the generalized Navier boundary condition (GNBC), and it is aimed 
at resolving the so-called contact line paradox of the movement of a two phase fluid on a solid wall. 
The reader is referred to, for instance, 48 , 49, 25 for a discussion of its derivation. Although there 
has been a lot of discussion and controversy around the validity of this boundary condition, see for 
instance [12 [61], we shall take the GNBC as a given and will not discuss its applicability and/or 
consequences here. 



2.2. Nondimensionalization. Here we present appropriate scalings so that we may write equa- 
tions (2.15|— (2.181 in non-dimensional form. Table 2.1 shows some typical values for the material 



ELECTROWETTING 



9 



Table 2.1. Physical Parameters at standard temperature (25° C) and pressure (1 
bar), taken from |35j. A Farad (F) is C 2 /J. For drinking water, K is 5 • 10~ 4 to 
5-Hr 2 . 



Parameter 


Value 


Surface Tension 7 


(air/water) 0.07199 J/m 2 


Dynamic Viscosity rj s 


(water) 8.68 • KT 4 , (air) 1.84 • 10" 5 Kg/m • s 


Density p a 


(water) 996.93, (air) 1.1839 Kg/m 3 


Length Scale (Channel Height) L s 


50 • lO" 6 to 100 • 10~ 6 m 


Velocity Scale U s 


0.001 to 0.05 m/s 


Voltage Scale V B 


10 to 50 Volts 


Permittivity of Vacuum e vac 


8.854- 10" 12 F/m 


Permittivity e s 


(water) 78.36 • e va c, (air) 1.0 • £ vac 


Charge (Regularization) Parameter A 


0.5 J-m 3 /C 2 


Mobility M s 


0.01 m 5 /(J-s) 


Phase Field Parameter a 


0.001 J-s/m 2 


Electrical Conductivity K s 


(deionized water) 5.5 • 10 -6 , 
(air) 0.0 C 2 /(J • m • s) = Amp/ (Volt • m) 



parameters appearing in the model. Consider the following scalings: 



p 


= p/ p s (choose p s ), 


V 


= 77/773 (choose 7?s), 


P = 




P s 


= Vs/L s , 


p 


= P/Ps, 


Pa 


= PsU 2 s , 


u = 


u/U s (choose U B ), 


X 


= x/L s (choose L s ), 


t 


= t/t s , 


is 


= L s /U s , 


P = 


P/Ps, 


Ps 


= l/L a , 


q 


= q/q s , 


q s 


= v s /x, 


V = 


V/Vs, (choose V s ), 


i 


= e/e a , 


5 


- S/L s , 


M 


- M/M s , 


K = 


K/K s , 


Ca 


VsU s 

1 


Re 


_ p s U s L s 

Vs 


We 


_ p s U*L s 

• 

7 


Bo = 


e s Vs 2 
L s ~f ' 


Je 


q s V s ' 


s P 


1 

a/t s ' 


M 


7 M S 
LIU* 


K = 


VsK s 
L s q s U s ' 




q s L 2 s 

V s es' 



where Ca is the capillary number, Re is the Reynolds number, We is the Weber number, Bo is 
the electro-wetting Bond number, 2e is the ratio of fluid forces to electrical forces, <Sp is the ratio 
of surface tension to "phase field forces," Mo is a (non-dimensional) mobility coefficient, Kq is a 
conductivity coefficient, and Ch is an electric charge coefficient. 
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Let us now make the change of variables. To simplify notation, we drop the tildes, and consider 
all variables and differential operators as non-dimensional. The fluid equations read: 

7 S(u)) + Vp = - — gV (V + q) + ip'(0)0 t u, in fi, 

in f2, 

on r, 

---^(Q'fM + ^d^, onT. 

The phase-field equations change to (again dropping the tilde) 

b t + u-V(j) = M o V-(M(0)V/i), inO, 



(D(pu 


) 1 


Dt 


Re 


Vu = 

< 


o, 


u-n = 


0, 




»?S(u) 



/i = Qw'(^) - SAA - Bo* £ '(4>)| W| 2 + Wei//(^)|u| 



m 



0. 



>< + u T d T <£ + S P {Q' f M + 5d n (j>) = 0, d n n = 0, on r. 

Performing the change of variables on the charge transport equation gives 

U + V-(qu) = K V {K{4>)V (q + V)) , in Q, 
jn-V (q + V) = 0, onT. 

Lastly, for the electrostatic equation we obtain 

-V- ((£)W) = C mX n, inQ*, 
V = V /V s , on d*n\ 

d n v = o, onao*nr. 

where £*(</>) has been normalized by e s . 

To alleviate the notation, for the rest of our discussion we will set all the nondimensional groups 
(Ca, Re, We, Bo, Ie, £p, Mo, Kq and Ch) to one. If needed, the dependence of the constants on 
all these parameters can be traced by following our arguments. Moreover, we must note that if a 
simplification of this model is desired, then these scalings must serve as a guide to decide which 
effects are dominant. 



2.3. Tangential Derivatives at the Boundary. As we can see from (2.151 and (2.16), our model 
incorporates tangential derivatives of the phase variable 4> at the boundary T. Unfortunately, in the 
analysis, we are not capable of dealing with these terms. Therefore, we propose some simplifications. 

The first possible simplification is simply to ignore the terms that contain this tangential deriva- 
tive; see [2D]. However, it is our feeling that the presence of them is important, specially in dealing 
with the contact angle in the GNBC. 

A second possibility would be to add an ad hoc term of the form Ar^ on the boundary condition 
for the phase variable, where by Ap we denote the Laplace-Beltrami operator on T. A similar 
approach has been followed, in a somewhat different context, for instance, by Priiss et al. [47] and 
Cherfils et al. [14] . However, this condition might lead to lack of conservation of <f>, which is an 
important feature of phase field models based on the Cahn Hilliard equation. 

Finally, the approach that we propose is to recall that, in principle, the phase field variable 
must be constant in the bulk of each one of the phases and so d T <fi ~ there. Moreover, in the 
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sharp interface limit this tangential derivative must be a Dirac measure supported on the interface. 
Therefore we define a function 



(2.20) 



1 i 



where 5 is non-dimensional, 



and replace all the instances of d T <j) by ip (</>)■ 



2.4. Contact Line Pinning. Simply put, the contact line pinning (hysteresis) is a frictional effect 
that occurs at the three-phase contact line, and is rather controversial. We refer the reader to 
[Ml R)3"] for an explanation about its origins and possible dependences. Let us here only mention 
that, macroscopically, the pinning force has a threshold value and, thus, it should depend on the 
stress at the contact line. It is important to take into account contact line pinning since, as observed 
in [5H , it is crucial for capturing the true time scales of the problem. 

We propose a phenomenological approach to deal with this effect. From the GNBC, 

f3u T + v S(u) nT = 7 (Q' f3 (ct>) + £cW) ^(0), 

we can see that, to recover no-slip conditions, one must set the slip coefficient j3 sufficiently large. 
On the contrary, when f3 is small, one obtains an approximation of full slip conditions. A simple 
dimensional argument then shows that f3 = rji, where £ has the dimensions of inverse length. 
Therefore, we propose the slip coefficient to have the following form 



= TlWifaS), 



where 



£(0,S) = 



6, 



W > 2' 

|0| < - and |S(u) nT | «T p 



|0| < -, and |S(u) r 



|0| < -, and |S(u) nT | »T p 



where 6 is the non-dimensional transition length. For the purposes of analysis, we face the same 



difficulties in this expression as in [2.3 Ergo, we will use this to model pinning in the numerical 
examples, but leave it out of the analysis. 



3. Formal Weak Formulation and Formal Energy Estimate 



In this section we obtain a weak formulation for problem ( 2.15 1— ( 2.18) and show a formal energy 



estimate, which serves as an a priori estimate and the basic relation on which our existence theory 
is based. 



3.1. Formal Weak Formulation. To obtain a weak formulation of the problem, we begin by 
multiplying the first equation of (2.15) by 0, the second by /2 and integrating in fi. After integration 



by parts, taking into account the boundary conditions, we arrive at 
(3.1a) (0 t , 0) + <u-V0, 0) + (M(0)V M , V0) = 0, 
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and 

(3.1b) (»,P) = Z (W'(ct>),p) + 7 <W,V/2) - \ (e'(ct>)\VV\ 2 ,p) + i <p'W0|u| 2 , /2> 



a [0 t + u r ip((p),fj] + 7 [e} s (0), a*] 



Multiply the first equation of (2.16) by w such that w-n|r = 0, the second by p and integrate in £1. 
Integration by parts on the first equation, in conjunction with the boundary conditions and (2.201, 
yields 

- (V-fafo)S(u)), w) = fa(#S(u), S(w)) - fofa)S(u) ni w x ] 

= (r,(0)S(u),S(w)) + [/3(0)u T ,w T ] - 7 [& f M + ^,w T ^(^)] 
= (ij(0)S(u), S(w)) + [/3(^)u T) w T ] + a [</h + w T ^(^)] , 



where we used the third equation of (2.15|. With these manipulations we obtain 



(3.2a) 



gW) ") 



,w) + (j 7 (0)S(u),S(w)) - (p,V-w) + [/3(0)u x ,w T ]+a[u x V(0),w T V(^)] 

= (^V0, w> - (qV(Aq + V), w> + - (p'{^)<j> t U, w)~a [&^($, W T ] , 



for all w, and 
(3.2b) 



(p,V-u> =0, 

for all p. Multiply (2.17) by r and integrate in f2 to get 

(3.3) (q t , r) - (<?u, Vr) + (if (^)V(Ag + V), Vr) = 0. 

Let be a function that equals zero on d*Sl* . Multiply the equation for the electric potential 
(2.18) by W, integrate in fl* to obtain 

(3.4) (s*(cf>)VV,VW) nir = (q,W). 

Given the way the model has been derived, it is clear that an energy estimate must exist. Before 
we obtain it let us show a comparison result a la Grdnwall. 

Lemma 3.5 (Gronwall). Let f,g,h,w : [0, T] — > K be measurable and positive functions such that 

(3.6) f(t) 2 + f g(s)ds<h(t) + f f(s)w(s)ds, Vtg[0,T]. 

Jo Jo 

Then ^ 

sup f{s)' 2 + ll .g(s)ds<4 sup h{s)+4T I w 2 (s)ds, Vie[0,T] 
«e[o,T] 2 J se[o,T] Jo 

Proof. Take, in (3.6), t = to, where 

to = argmax{/(s) : s € [0,T]}, 

then 

/(*o) 2 + / g(s) ds < max h(s) + f(t Q ) [ ° to(a) ds < max fc(s) + \f{t a ) 2 + 



s£[0,T] 



se[o : T] 



w(s) ds 



Canceling the common factors, applying the Cauchy Schwarz inequality on the right and taking the 
supremum on the left hand side we obtain the result. □ 
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Remark 3.7 (Exponential in time estimates). The main advantage of using Lemma 3.5 to obtain 
a priori estimates, as opposed to a standard argument invoking Gronwall's inequality, is that we 
can avoid exponential dependence on the final time T. 

The following result provides the formal energy estimate. 



Theorem 3.8 (Stability). If there is a solution to ( 2.15 ) — ( 2.18 1 7 then it must satisfy the following 
estimate 



(3.9) 



sup 

«e(o,T] Un 

T 



P (0)|u| 2 + V+ 7 ( ^|V0| 2 + 9W(0) 



+ / j£ *(0)|vvr + 7 / QfM 



UCl 



[t7(0)|S(u)| 2 + M(4>)\Vp\ 2 + K(<t>)\V(\q + V)\ 2 } + / [/3(</>)|u.| 2 + a|0 t + uW#)| 



< 



P (0)|u| 2 +q 2 + 7 -|V0| 2 + -W(0) + -\V Q \ 



+7 / Ofs(4>) 



+ sup 

4=0 se[o,r] Un 



e M \VV \ 



cT 



£*(0)|W| 2 +£ M |W o | 2 ) 
1 



A 



Vo\ z (t) 



A 



|Vn t | 2 



where c does not depend on T. 

Proof. We first deal with the Navier Stokes and Cahn Hilliard equations in a way very similar to 
Theorem 3.1 of [S3]. Set w = u in (3.2a| and notice that 

1 d 



D(pu) 
Dt 



2 dt 



P|u| 2 , 



because 

We obtain 

d 1 



(3.10) 



dt 2 



D (P U ) , S r, s 

-g^ = a(au) t + pu-Vu + -V(pu)u. 



p\u\ 2 + / 77|S(u)| 2 + / /3(0)|u T | 2 + a / |u T V(0)| 2 = W,u> 



(qV(\q + V),u) + - (p' (<!>)&, |u| 2 ) - a [cf> t u T , ^{cj>)] 



Set 4> = M m (3. la | to get 
(3.11) 



(/*,&> + W,u)+ / M(^»)|V/i| 2 = 



Set p - 
(3.12) 



^ in (3.1b ) to write 
d 



. M) = -7 



df 



J2 



+ -< £ '(^ t ,|vy| 2 ) 



i 



<p'(0)0 t ,|u| 2 )-a / (^) 2 -a[^,u^(^ 
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Add (3.10), (3.11) and (3.12) to arrive at 



(3.13) 



df 



^)|S(u)| 2 



+ / p{4>)\u T \ 2 + M(0)|V/i| 2 +a / (0 t + u^(0)) 2 = -(gV(Ag + y),u) + -( e '(0)0 t ,|VT/| 2 ). 



We next deal with the electrostatic equations. Set r = Xq + V in (3.3) to get 
A d 



(3.14) 



2 dt 



q 2 + (V,q t )~(qV(Xq + V),u)+ / X(0)| V(Ag + V)| 2 = 



Take the time derivative of (3.4) and set W = V — Vo, where by Vo we mean an extension of Vo to 
O*. We obtain 



(3.15) / dt(e* (</>)) \VVf 



1 



e*{<t>)d t (|VV| 2 ) = V) - (q u Vo) + ($(e*(0) W), VV&) n . • 



Add (3.13), (3.14) and (3.15) and recall that £*(</>) is constant on fi* \ £1 We thus obtain 



d 

dt 



i /9 (0)|u| 2 + ^ 2 + 7 ^|V0| 2 + ~w(0 



[t7(0)|S(u)| 2 + M(0)|V M | 2 + ^(0)|V(Ag + V)| 2 ] + y [/3(0)K| 2 + + u T ^)\ 2 ] 

= (dt(e*((f>)VV), VV^o> n *- (ft, Vb) • 



Integrate in time over [0, t], with < t < T and integrate by parts the right hand side. Repeated 
applications of the Cauchy-Schwarz inequality give us 



^)|u| 2 + ^ 2 + 7 ^|V^| 2 + ^WW 



1 



< 



+ y^ -^(0)| W | 2 + 7 y e /a (^) 

[t7(0)|S(u)| 2 + M(0)|Vm| 2 + ^(0)|V(Ag + V)| 2 ] + / / [/3(<£)|u T | 2 + a|& + u T V#)| 2 ] 



o -T 



^(0)|u| 2 + ? 2 + 7 ^|V^| 2 + iwW 



+ / £ *(0)|vv-| 2 + 7 / e /a (0) 

Jn* Jr 



t=o 



e M (\VV \ 2 (t) + |VVb| 2 (0)) + 
si* Jn 



emIvv^i 2 + ^ / |y , t ' 2 



^|Vo| 2 W + ^|Vo| 2 (0) 



o Un 



1/2 



A 



1 



^)|Vlf 



1/2 



where em is the maximal value of the function £*(</>). 
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Finally, if we set 
1 



/(*) = 



2 ^)|u| 2 + ^ 2 + 7 (j|V0| 2 + - s W(cf>) 



+ J ^)|W| 2 + 7 ^e /s (<^(/). 



g(t)= / [r;(0)|S( U )| 2 + M(0)|V/z| 2 + ^)|V(A 9 + y)| 2 ] + / [/3(0)|u T | 2 + a\cf> t + u T ^)\ 2 } (i), 
Jn Jr 



h(t) 



w(t) 



^p(0)|u| 2 + g 2 + 7 (||V0| 2 + i>V(0) 



em (|Vt^ | 2 (t) + |W | 2 (0)) - 
n* Jn 

1/2 



+ / £ ^)|w| 2 + 7 / e /s (0) 

Jn* Jr 
1 



X |Vo| 2 W+ 2 |Vo| 2 (0) 



e M NV A z + 



\Vo,, 



then an application of Lemma 3.5 gives the desired estimate 



□ 



4. The Fully Discrete Problem and Its Analysis 
In this section we introduce a space-time discrete problem that is used to approximate the 



electrowetting problem (3.1 H3.4). Using this discrete problem, and the result of Theorem 3.8 we 



will prove that a time-discrete version of our problem always has a solution. Moreover, in Section[5] 
we will base our numerical experiments on a variant of the problem defined here. 



4.1. Definition of the Fully Discrete Problem. To discretize in time, as discussed in §!.![ we 
divide the time interval [0, T] into subintervals of length At > 0. Recall that the time increment 



operator was introduced in (1.4) and the time average operator (•) in ( |l.5[ ). 

To discretize in space, we introduce a parameter h > and let W/j C Hl(Q*), Qh C 1 (il) 
X h C V and M h 

condition between the spaces and 



C L? = (Q) be finite dimensional subspaces. We require the following compatibility 



(4.1) 



W h \n G 



VW h G 



Moreover, we require that the pair of spaces (X^M^) satisfies the so-called LBB condition (see 
HH HI] ) > that is, there exists a constant c independent of h such that 



(4.2) 



c\\ph\\L 2 < SUp 



Vph G 



Finally, we assume that if Y is any of the continuous spaces and Y^ the corresponding subspace, 
then hi < hi implies Y/, 2 C Y/^. Moreover, the family of spaces {Yh}h>o, is "dense in the limit." 
In other words, for every h > there is a continuous operator Z/j : Y — > Y^ such that when h — > 

\\y-IhVh->0, VyeY. 

The space will be used to approximate the voltage; Qh the charge, phase field and chemical 
potential; and X^, the velocity and pressure, respectively. Finally, to account for the boundary 
conditions on the voltage, we denote 



k + l\ 



v k+1 - 
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Remark 4.3 (Finite elements). The introduced spaces can be easily constructed using, for instance, 
finite elements, see [TTJ [3TJ [TB] for details. The compatibility condition (4.1) can be easily 
attained. For instance, one can require that the mesh is constructed in such a way that for all cells 
/C in the triangulation Tk, 

K n Q ^ & k n (n* \ n) = 0, 

and the polynomial degree of the space Qh is no less than that of W^. Finally, we remark that the 
nestedness assumption is done merely for convenience. 

The fully discrete problem searches for 

{V hA t ~ V ,AU<lhAt,<f>hAU HhAt,UhAuPhAt} C W h x(jjx X h X M h , 

that solve: 

Initialization: For n — 0, let q^, 0° and u° be suitable approximations of the initial charge, 

phase field and velocity, respectively. 
Time Marching: For 0<n<iV— lwc compute 

(v ^+l )g „+l^n+l )M n+l )U „ + l )p „+l ) e Wh (yn+^ x Q3 x Xft x Mfc> 

that solve: 

(e*(<f> n h +1 )W£ + \VW h ) nt = (q h l+1 , W h ) , VW h e W h , 



(4.4) 
(4.5) 

(4.6) 



At 



r h } - (« + \ Vr h ) + <Jf(0*)V (A^ l+1 + , Vr„) = 0, Vr* e Q h , 



k n+l 



At 



$ h ) + « +1 -v^,4) + (M(<^)V< + \ vfa) = o 



&fc G Qh 



(4.7) = j (W'M)+^ +1 ,^)+7*(V^ +1 ,VM h ) 



■ 7 [&' fs m + B ^h +1 ^h] Vfih e Q h , 



where we introduced 
(4.8) £{<Pi,<P2) 



e' (scp! + (1 - s)(p 2 ) ds, 



(4.9a) 



At 



n(th n )n n \ 1 



< ?7 (^)S«+ 1 ),S(w, l )) - <K +1 , Vw,) + [/3(^X+ 1 ) w ftT ] +a [u£V(^),w AT ^)] 



— a 



At 



Vw,, e 
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(4.9b) 



Remark 4.10 (Stabilization parameters). Notice that, in (4.7 1, we have introduced two stabiliza- 
tion parameters, namely A and B. Their purpose is two-fold. First, they will allow us to treat 
the nonlinear terms explicitly while still being able to mantain stability of the scheme, see Propo- 
sition [47T4] below. Second, when studying convergence of this problem, the presence of these terms 
will allow us to obtain further a priori estimates on discrete solutions which, in turn, will help in 



passing to the limit, see Theorem 6.7 We must mention that, this way of writing nonlinearities 
is related to the splitting of the energy into a convex and concave part proposed in [65] . See also 



Remark 4.11 (Derivative of the permittivity). Notice that (4.8 1, i.e., the definition of the term £ , 
is a highly nonlinear function of its arguments (unless e is of a very specific type). As the reader 
has seen in the derivation of the energy law (Theorem 3.8 ), the treatment of the term involving the 
derivative of the permittivity is subtle. In the fully discrete setting this is additionally complicated 
by the fact that we need to deal with quantities at different time layers. The reason to write the 
derivative of the permittivity in this form is that 

f gOgi) -e(v?2) 

ly(¥»i)i <pi = <P2, 

which will allow us to obtain the desired cancellations. 



The following subsections will be devoted to the analysis of problem (4.4)— (4.9). For convenience, 
we define 

Tirh n 

^:=^+u^r 1 ). 



4.2. A Priori Estimates and Existence. Let us show that, if problem ( 4.4 H 4.9 1 has a solution, 
it satisfies a discrete energy inequality similar to the one stated in Theorem |3.8[ To do this, we 
first require the following formula, whose proof is straightforward. 

Lemma 4.12 (Summation by parts). Let {/"}^tTo an ^ {d n }ri=o be sequences and assume f^ 1 = 
g~ x = 0. Then we have 

(4.13) 



m— 1 

E 

n=0 



m— 1 „m — 1 



{dg n )f n = f m-, g 



m-2 

E 

71 = 



S"(f/ n+i )- 



Proposition 4.14 (Discrete stability). Assume that the stabilization parameters A and B are 
chosen so that 



(4.15) 



A>lsn V W"(0, B>isupe?,(0. 



The solution to ( 4.4 1— (4.9), if it exists, satisfies the following a priori estimate 



(4.16) ||u^ A t||£^(L 2 ) + II^U/iAtllljV^L 2 ) + ll u /iAt||£2( V ) + || QhAt (L 2 ) + II ^QhAt II f,V 2 (L 2 ) 
+ ||V0( l Aj£^(L 2 ) + l|VO<KAjl,i/ 2 (L 2 ) + l|W(</>hAi)^^ 

+ II^At|l<? 2 (L 2 (r)) + ||©/ s (</>ftAt)ll^(Li(r)) + HVM/iAtll^fL 2 ) + II V(Ag + "l / ), lA t||£ 2 (L 2 ) < c, 
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where we have set pP h = 0, Vfi = for convenience in writing (4.161. The constant c depends on 
the constants 7, S, a, the data of the problem u°, ^>°, V"o,At and T, but it does not depend on 
the discretization parameters h or At, nor the solution of the problem. 



Proof. We repeat the steps used to prove Theorem 



i.e., set w;,. 



2Atu™ +1 in (|4.9a 

jn+1 



in ( |4.9b| ), cj> h = 2At^ +L in ( |4.6D , = -2u<^ +i in ((^Tj) and r h = 2At(A<^ +i + V^* 1 ) in ([4^Sj> . To 
treat the time-derivative terms in the discrete momentum equation, we use the identity 



, Ph 



Pl +1 



2< + 



^ +1 )K +1 | 2 



P {^ h )\K\ z + p{^K 



n+l\2. 



see [311132) . To obtain control on the explicit terms involving the derivatives of the Ginzburg-Landau 
potential W and the surface energy density notice that, for instance, 



for some £. Choosing the stabilization constant according to (4.15) (cf. 59, 58, 60, 53]), we deduce 
that 

' (W'{^ h )+A^l +1 )^l +1 > 



Adding (|4.5[)-(|4.9[) yields, 
(4.17) o|| ( t(^+ 1 )<+ 1 ||^ 



11- ilj^^ 



2 A 

6 



n - in 2 2 



n+l||2 



n+l 



L 2 T) 



27 



06 - 



0W(^ +1 ). 

! +ii^r i ii 2 L 2 )+7^(9iiv0r i n L2 



r 

n+l 



2 

L 2 



2At 



n+l- 



L 2 



^(«)V(A< 



n+l 



n+l\ 



2 

L 2 



A«+l 



n+l 



At 



L 2 (r) 



Take the difference of ( 4.4 ) at time-indices n + l and n to obtain 



and set ^ = 2(V h n+1 - ^ " +i ). In view of Jl.6|) we have 



20 (£*(C +1 )VV™ +1 ) -W h n+1 = (e*(0" +1 )|VV^ +1 l 2 ) 

+ e*(^)|V0K 1+1 | 2 +5( £ *(C +1 ))|V^r 

whence 

2 



(4.18) 



L 2 (a* 



n+l 
h 



2 




H 

L 2 (n*) 


-/ 



o £ *(^ +1 )|vy, n+1 l 2 



2 - 2 v^ 1 ) + 2 (0 ( £ *(0r 1 )vv^ +1 ) , w n+1 ) n . . 



Add (4.17) and (4.18). Notice that, since the permittivity is assumed constant on fi*\f2, on the 
left hand side of the resulting inequality we have the following term 



in+l\ 



)-£(<Pt + \<t>i)wr l )\vv£+ i \ 2 =o 
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where we used the definition of E, see (4.8 1 and Remark 4.11 Therefore, we obtain 

(4.19) o|K<^ +1 k h 



kn+l\,,n+l 1 1 2 

II 2 



n+1 1|2 ^ 



A ! f>IK +1 " 2 



+ j6 (B||V^ 



n+l |2 
L 2 



L 2 

2" 



2ll^ +1 Hi 2 



fW(C +1 ) 



+ 5 



m+l- 



71+1 



L 2 (fi*) 



n+1 



L 2 (n*) 



+ 27^ D6 /s 



+ 2At 



n+l* 



L 2 



M(^)V< 



n+l 1 1 2 
L 2 



^(K)v(A<+ 1 + y / i 



ra+l\ 



a 



L 2 



L 2 (r) 



m+l 



At 



< 



L 2 (r) 

rn+l\ 



-2<c)< +1 ,V "+ 1 ) + 2(c. ( £ *(^ +1 )VC +1 ),VKr +1 > n4 



Summing (4.19) for n — 0,...,m — 1, using summation by parts ( 4.13| ) (set yP h = 0, V® = 0), 
applying the Cauchy-Schwarz and weighted Young's inequality, we obtain the result. □ 



Remark 4.20 (Compatibility). Notice that condition (4.1) is needed to obtain the stability esti- 
mate, otherwise 2At(\q^ +1 + V£ + ) would not be an admissible test function for (4.5). 



The a priori estimate (4.16) allows us to conclude that, for all h > and At > 0, problem 



(4.4)-(4.9l has a solution. 



Theorem 4.21 (Existence). Assume that the discrete spaces satisfy assumptions (4.1) and (4.2), 
the stabilization parameters A, B are chosen as in Proposition \4-14\ Then, for all h > and 
At > 0, problem (4.4) -(4.9) has a solution. Moreover, any solution satisfies estimate (4.16). 



Proof. The idea of the proof is to use the "method of a priori estimates" at each time step. In 
other words, for each time step we define a map C n+1 in such a way that a fixed point of £ n+1 , if it 
exists, is a solution of our problem. Then, with the aid of the previously shown a priori estimates 
we show that C n+1 does indeed have a fixed point. 

We proceed by induction in the discrete time and assume that we have shown that the problem 
has a solution up to n. For each n = 0, N — 1, we define 

C n+1 : W h (V™ +1 ) xQlxX h x M h -+ W h (V£ +1 ) xQ^x X h X M h , 



{V h ,qh,<t>h,Hh,u- hl p h ) 

where the quantities with hats solve 
(4.22) 



(4.23) 



Qh - K 
At 



£*Wv%, VW V^ = <®" Wh) ' yWh e Wh ' 



, r h \ - {q h u h , Vr h ) + (if(^)V (a% + V^) , Vr fc ) = 0, Vr/, e Q h , 



(4.24) 



At 



, 4 ) + (iifc-VC ^) + (M(^)VAa, V&) = 0, 



5>h € 
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(4.25) (p h ,p h ) 



| (W'(^) +A(<p h - ,fl h ) + 7 ^(v4, V/2 h ) + - (p'(<Pl)u%-u h ,(i h ) 



^(£{4> hl 4>l)WV h -WV h ,JXh) + a 



At 



+ 7 [& fa {ft>D +B(cj> h - 4>l) , p h ] V/ifc G Q h) 



(4.26) 



(p(^)<-Vu ft ,wO 



1 



:v-(p(«w)fi„,w h ) 



At ' "'7 ' • ""' ' 2 



— U/^W/, 



At 



Vwh G 



(4.27) 



(p h ,Vu h )=0, Vp h £M h . 



Notice that a fixed point of C n+1 is precisely a solution of the discrete problem (4.4 )— (4.9 1. 
To show the existence of a fixed point we must prove that: 

• The operator C n+1 is well defined. 

• If there is a X = (Vh, qh, §K)P>hi Uh,Ph) for which X = uj£ n+1 X, for some uj G [0, 1], then 
(4.28) \\X\\<M, 

where M > does not depend on ^ or u. 

Then, an application of the Leray-Schauder theorem [22j [66] will allow us to co nclud e. Moreover, 
since a fixed point of C n+1 is precisely a solution of our problem, Proposition 4.14 gives us the 
desired stability estimate for this solution. 

Let us then proceed to show these two points: 
The operator C n+1 is well defined: Clearly, for any given <\>h, and qh, the system (4.22 )— ( 4.23 1 is 



positive definite and, thus, there are unique Vh and qh- Having computed Vh and qh we then 



notice that (4.26) and (4.27) are nothing but a discrete version of a generalized Stokes problem. 



Assumption (4.2) then shows that there is a unique pair (iih,Ph)- To conclude, use (Vh, q^, Uh,Ph) 
as data in (4.24) and (4.25). The fact that this linear system has a unique solution can then be 



seen, for instance, by noticing that the system matrix is positive definite. 

Bounds on the operator: Notice, first of all, that one of the assumptions of the Leray-Schauder 
theorem is the compactness of the operator for which we are looking for a fixed point. However, 
this is trivial since the spaces we are working on are finite dimensional. Let us now show the 
bounds noticing that, at this stage, we do not need to obtain bounds that are independent of h, 
At or the solution at the previous step. This will be a consequence of Proposition |4.14| Let us 
then assume that for some X — (Vh,qh,<f>hi Phi u hiPh) we have X = ui£ n+1 X. Notice, first of all, 
that if to = then X = and the bound is trivial. If oj e (0, 1], the existence of such element can 
be identified with replacing, in ( |4.22| )- (j4.27 |, (VhAhAh,fah,&h,Ph) b y u>~ \V h , qh, <t>h,Ph, U- h,Ph) - 
Having done that, set w h = 2A tu h in ( |4.26[ ), r h = 2At(\q h + V h ) in ( |4.23[ ), <j> h = 2Atp h in \A.2A\ 
and ph = 2{4>h — 4>h) in (4.25). Next we observe that, by induction, the equation has a solution 



at the previous time step, therefore there are functions that satisfy (4.4) for time n. Multiply this 
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identity by u> and subtract it from (4.22). Arguing as in the proof of Proposition 4.14 we see that 



condition (4.15) implies that to obtain the desired bound we must prove estimates for the terms 



(p'^Ku,,^), (n h ,<%), (qlV h ) 



which are, in a sense, the price we are paying for not being fully impicit. All these terms are linear 
X and, thus, can be easily bounded by taking into account that we are in finite dimensions and 
that the estimates need not be uniform in h and At. □ 



5. Numerical Experiments 

In this section we present a series of numerical examples aimed at showing the capabilities of 
the model we have proposed and analyzed. The implementation of all the numerical experiments 
has been carried out with the help of the deal . II library [71 [5] and the details will be presented in 



Let us briefly describe the discretization technique. Its starting point is problem (4.4)-(4.9) 
which, being a nonlinear problem, we linearize with time-lagging of the variables. Moreover, for 
the Cahn Hilliard Navier Stokes part we employ the fractional time-stepping technique developed 
in |53j . In other words, at each time step we know 

(F^C/M,PM G W fe (F ") x Ql x X, x M 2 , 

with £jj := and, to advance in time, solve the following sequence of discrete and linear problems: 

• Step 1: 

Potential: Find V," +1 £ Wfc(V n+1 ) that solves: 

(e*(4%)VV£ +1 , VW h ) nt = (ql,W h ) , VW h £ W h , 
Charge: Find q£ +1 £ Q h that solves: 

• Step 2: 

Phase Field and Potential: Find £ Q h that solve: 

4> h \ + (u£-VC 4) + <M(^)V^ +1 , V^) = 0, £ Q h , 



in+1 



At 



■ 7 [ej 



fsWh) 



Step 3: 
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Velocity: Define p\ = p% + then find u£ +1 G X h such that 

^rv^-^K w ^ + (pW)u „ Vu „ +1; wh) + 1 (vw « KK+ i )W/i) 

+ (r / (^)S(u;;+ 1 ),S(w h ))-(p»,Vw, l )+ [«)u^\w^] +a[u»+V(^),w,^(^)] 

1 / ,, 



1 / t)(A" +i \ 



-\ i n-\-l 

Step 4: 

Penalization and Pressure: Finally, and J>^ +1 are computed via 



where p := min{ j oi,p2} and 



Remark 5.1 (CFL). A variant of the subscheme used to solve for the Cahn Hilliard Navier Stokes 
part of our problem was proposed in [53| and shown to be unconditionally stable. In that refer- 
ence, however, the equations for the phase field and velocity are coupled via terms of the form 
(u^ +1 -V<^, 4>h)- If we adopt this approach, coupling steps 2 and 3, and assume that the permit- 
tivity does not depend on the phase, it seems possible to show that this variant of the scheme 
described above is stable under, 

At < cSh. 

On the other hand, if we work with full time-lagging of the variables, then it is possible to show 
that the scheme is stable under the, quite restrictive, assumption that 

At < c6 2 h 2 . 

To assess how extreme these conditions are one must remember that, in practice, it is necessary to 
set h — 0{8). Nevertheless, computations show that these conditions are suboptimal and just a 
standard CFL condition is necessary to guarantee stability of the scheme. 

5.1. Movement of a Droplet. The first example aims at showing that, indeed, electric actuation 
can be used to manipulate a two-fluid system. The fluid occupies the domain — (—5, 5) x (0, 1) 
and above and below there are dielectric plates of thickness 1/2, so that il* = (—5, 5) x (—1/2, 3/2). 
A droplet of a heavier fluid shaped like half a circle of radius 1/2 is centered at the origin and 
initially at rest. To the right half of lower plate we apply a voltage, so that 

V = VooXD, D = \{x, y) G R 2 : x > 0, y = -} 



The density ratio between the two fluids is p\j p-2 — 100, the viscosity ratio 771 / 772 = 10 and the 
surface tension coefficient is 7 = 50. The conductivity ratio is K1/K2 = 10 and the permittivity 
ratio ei/ 82 — 5 and Sd/s2 — 100. We have set the mobility parameter to be constant M = 10 -2 , 
and a = 10~ 3 . The slip coefficient is taken constant j3 = 10, and the equilibrium contact angle 
between the two fluids is 6 S — 120°. The interface thickness is S = 5 ■ 10~ 2 and the regularization 
parameter A = 0.5. The applied voltage is Vqq = 20. 
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The time-step is set constant and At = 10 -3 . The initial mesh consists of 5364 cells with two 
different levels of refinement. Away from the two-fluid interface the local mesh size is about 0.125 
and, near the interface, the local mesh size is about 0.03125. As required in deal. II, the degree of 
nonconformity of the mesh is restricted to 1 i.e., there is only one hanging node per face. Every 10 
time-steps the mesh is coarsened and refined using as, heuristic, refinement indicator the L 2 -norm 
of the gradient of the phase field variable (f>. The number of coarsened and refined cells is such that 
we try to keep the number of cells constant. 

The discrete spaces are constructed with finite elements with equal polynomial degree in each 
coordinate direction and 

degW h = l, degQ ft = 2, degX fe = 2, dcgM h = 1, 

that is the lowest order quadrilateral Taylor-Hood element. No stabilization is added to the mo- 
mentum conservation equation, nor the convection diffusion equation used to define the charge 
density. 



Figure 5.1 shows the evolution of the interface. Notice that, other than adapting the mesh so 
as to resolve the interfacial layer, no other special techniques are applied to obtain these results. 
As expected, the applied voltage creates a local modification variation in the value of the surface 
tension between the two fluids, which in turn generates a forcing term that drives the droplet. 

5.2. Splitting of a Droplet. One of the main arguments in favor of diffuse interface models is 
their ability to handle topological changes automatically. The purpose of this numerical simulation 
is to illustrate this by showing that, using electrowetting, one can split a droplet and, thus, control 
fluids. Initially a drop of heavier material occupies 



The material parameters are the same as in |5.1| To be able to split the droplet, the externally 
applied voltage is 

D=[(x,y)eR 2 :\x\>^, V = ~\ 

Figure [5~2| shows the evolution of the system. Notice that, other than adapting the mesh so as 
to resolve the interfacial layer, nothing else is done and the topological change is handled without 
the necessity to detect it or to adapt the time-step. 

5.3. Merging of Two Droplets. To finalize let us show an example illustrating the merging of 
two droplets of the same material via electric actuation. The geometrical configuration is the same 
as in §5.2| In this case, however, there are initially two droplets of heavier material, each one of 
radius 0.5 and centered at (—0.7,0) and (0.7,0), respectively. The material parameters are the 



same as in £5.2 except the interfacial thickness, which is set to 5 = 10 . We apply an external 
voltage so that 

D = [(x,y)eR 2 : \x\<±, V = -\ 

To be able to capture the fine interfacial dynamics that merging possesses, we set the initial 
level of refinement to 4, with 3 extra refinements near the interface, so that the number of cells is 
48, 696 with a local mesh size away of the interface of about 0.02875 and near the interface of about 
6 • 10~ 3 . This amounts to a total of 147,249 degrees of freedom. The time-step, again, is set to 
At = 1Q- 3 . 
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t=Q*lQA-3 ' t=200*iOA-3 



t=400*i0A-3 Z t=600*i0A-3 



t=800*l6A-3 t=1000'1QA-3 




t=1200*10A-3 t=1400*10A-3 



Figure 5.1. Movement of a droplet under the action of an external voltage. The 
material parameters are p\j pi = 100, 771/772 = 10, 7 = 50, K1/K2 = 10, £i/e2 = 5, 
e D /e 2 = 100, M = 10~ 2 , a = 10™ 3 , /3 = 10, 9 S = 120°, 5 = 5- 10~ 2 , A = 0.5 and 
V 00 = 20. The interface is shown at times 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2 and 1.4. 
Colored lines are used to represent the iso-values of the voltage. The black dotted 
line is the position of the interface at the beginning of the computations. 
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f =0* 1 6" A -3 t^250']QA-3 t=5ttHQA-3 




t=310Q'10A-3 t=3250"10"-3 t=350Q'10A-3 



Figure 5.2. Splitting of a droplet under the action of an external voltage. The 
material parameters are px/ p% = 100, r\\jr\i = 10, 7 = 50, KxjK^ = 10, £1/^2 = 5, 
e D /e 2 = 100, M = 10~ 2 , a = 10" 3 , /3 = 10, 6 S = 120°, 5 = 5- 10" 2 , A = 0.5 and 
Vqq = 20. The interface is shown at times 0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 
2.0, 2.25, 2.5, 2.75, 3.0, 3.02, 3.05, 3.10, 3.25 and 3.5. Colored lines are used to 
represent the iso- values of the voltage. The black dotted line is the position of the 
interface at the beginning of the computations. 

Figure [573] shows the evolution of the two droplets under the action of the voltage. Again, other 
than properly resolving the interfacial layer, we did not need to do anything special to handle the 
topological change. 

6. The Semi-Discrete Problem 

In §4.2| we showed that the fully discrete problem always has a solution and that, moreover, this 
solution satisfies certain a priori estimates. Our purpose here is to pass to the limit for h — > so 
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t=0*lQA-3 




t=2006' ; iQA-3 




t=5000*10A-3 "° " t=5500*10A-3 



Figure 5.3. Merging of two droplets under the action of an externally applied 
voltage. The material parameters are p\j pi = 100, 771/772 = 10, 7 = 50, K1/K2 = 
10, e x /e 2 = 5, e D /e 2 = 100, M = 10^ 2 , a = 10" 3 , P = 10, S = 120°, 5 = 10" 2 , 
A = 0.5 and Voo = 20. The interface is shown at times 0, 1, 2, 3, 3.3, 3.4, 3.5, 4, 
5 and 5.5. Colored lines are used to represent the iso- values of the voltage. The 
black dotted line is the position of the interface at the beginning of the computa- 
tions. 
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as to show that a semi-discrete (that is continuous in space and discrete in time) version of our 
electrowetting model always has a solution. 

Let us begin by defining the semi-discrete problem. Given initial data and an external voltage, 
we find: 

{V A t - Vb,At,8A*,0At,A»At,u A t,PAt} C Hl(Sl*) x H 1 ^) 3 x V x Lj =0 ( O ) 

that solve: 

Initialization: For n — 0, let q°, <f>° and u° equal the initial charge, phase field and velocity, 
respectively. 

Time Marching: For 0<n<iV— lwc compute 



G + K +1 x H^Slf x V x Lf =Q (O), 



(y n+i ,g" 

that solve: 

(6.1) (e*(0 n+1 )W n+ \W)^ = (g n+1 ,W), WGi^O*), 

(6.2) ^^p,r^-(g n u n+1 ,Vr) + (X(</.")V(Ag n+1 + y n+1 ),Vr) = 0, Vr G i/ 1 (Q) 

(6.3) ^^^,^ + ( u "+ 1 .V^,0) + (M(^)V^ +1 ,V^)=O, V^Gff 1 ^) 

(6.4) < M n+1 ,/2) = |(W'(0")+^" + \A2)+75<V0" + \VM) + i(p'(0")u".u" +1 ,ri) 



(£(0 n+1 ,0 n )|W" +1 | 2 ,/i} + , 



At 



(6.5a) / P^K^, ^"K w ^ + ^ p (^ )u «. Vu «+l + IvWflu-lu-^w 

+ ( f7 (0 n )S(u B+1 ),S(w)) - (p" +1 ,V-w) + [/3(^)< +1 ,w x ] + a ),w^(^)] 

1 / 7ifh n+1 
= (M n+1 Vf,w) - (q n V(\q n+1 + V n+1 ),w) + - (p'(cf> n )- 



At 



-u , w 



At 



, Vw G V, 



(6.5b) <p,Vu n+1 ) =0, VpeL^Sl). 

Remark 6.6 (Permittivity). Notice that, in our definition of solution, the test function for equation 



(6.4) needs to be bounded. This is necessary to make sense of the term 

<£(0"+\0")|W" +1 | 2 ,m), 

since £ is bounded by construction and V n+1 G i? 1 (f2). The authors of [2D] used a similar choice 
of test functions and showed, using different techniques, existence of a solution for their model of 
electrowetting in the case when the permittivity is phase-dependent. 
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Since the solution to the fully discrete problem ( 4.4 H 4.9 1 exists for all values of h > and 



satisfies uniform bounds, one expects the sequence of discrete solutions to converge, in some topol- 



ogy, and that the limit is a solution of problem ( 6.1 1— ( 6.5 1. The following result shows that this is 
indeed the case. 



Theorem 6.7 (Existence and stability). For all At > 0, problem ( 6.1 ) — ( 6.5 ( has a solution. 
Moreover, this solution satisfies an energy estimate, analogous to (4.16), where the constant c 



might depend on At and the data of the problem, but not on the solution. 



Proof. Theorem 4.21 shows the existence, for every h > 0, of a solution to the fully discrete problem 



(4.4)— (4.9 1 which, moreover, satisfies estimate (4.16). This estimate implies that, for every n, as 
h 0: 

• W (</)%) remains bounded in Since the modified Ginzburg-Landau potential is a 
quadratic function of its argument, this implies that there is a subsequence, labeled again 
4>h, that converges weakly in L 2 (fl). 

• remains bounded in L 2 . This, together with the previous observation, gives us a 
subsequence that converges weakly in and strongly in L 2 (f2). 

• The strong L? -convergence of (\>\ implies that the convergence is almost everywhere and, 
since all the material functions are assumed continuous, the coefficients converge also almost 
everywhere. 

• There is a subsequence of uj^ +1 that converges weakly in V and strongly in L 2 (il). 

• A subsequence of V£ — V$ converges weakly in il*(f2*) and hence strongly in L 2 (Q,*). 

• There is a subsequence of q^ +1 that converges weakly in L 2 (fl). Moreover, we know that 
K ((f>^)V (\q^ +1 + V? ) converges weakly. By the a.e. convergence of the coefficients and 
the L 2 -weak convergence of we conclude that V<7^ +1 must converge weakly and, 
thus, the convergence is weak in and strong in L 2 (f2). 

• The quantity </>^ +1 = -^j± t V U/^V-K'© remains bounded in L 2 (T), which implies that 

there is a subsequence of 0£ + that converges weakly in L 2 (T). 



V/iJJ remains bounded in L (fi). Moreover, setting /i;, = 1 in (4.7) and the observations 
given above, imply 



K + M)|< J(w'(^) + ^ +1 ,i) 



(p'(^K,< +1 ) + ( 



ifc+i 



which shows that J n remains bounded and, thus, /i^ remains bounded in i/ 1 (f2) and 
so there is a subsequence that converges weakly in and strongly in L 2 (£l). 



Finally, we use the compatibility condition (4.2 ) and the discrete momentum equation (4.9a 
to obtain an estimate on the pressure p r ^ + , 



c\\pI +1 \\lz < ^p(K)h™\\^ 



|v^|| L2 |K|| H1 |K 

U 2 (r) - 



Ih 1 - 



•«||^)ll^(r)||^ +1 ||L 2( r) + ||^ +l1 



u 



ra+li 



|L 2 



is(ur 



H 



L 2 



\1h \\h 



P(K)\\L 

l 2 + II«)IIl 



n+1 1 



Hi 



L 2 



+ ^l|p'(^)l|L~||O0r 1 ||L 2 ||<||v< 



At' 
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which, for a fixed and positive At, implies the existence of a L 2 -weakly convergent subse- 
quence. 

Let us denote the limit by 

{V A t ~ V o ,Au<lAu0AuVAuU AuPAt } c Hl(Sl*) x H^ilf x V x Lf =0 (Q). 
It remains to show that this limit is a solution of (|6.1|)-(|6.5|): 



(6.8) 



Equation (6.1 ): Notice that if we show that, as h — > 0, the sequence V^ +l converges to V n+1 
strongly in Hl(Cl*), then the a.e. convergence of the coefficients implies 

(e*(^ +1 )V^ l+1 , W) , -> (e*(0 n+1 )W B+1 ,VW> nik . 

Let us then show the strong convergence by an argument similar to that of |20l pp. 2778]. 
For any function V £ we introduce the elliptic projection VhV € Wh(V) as the 

solution to 

(WV h V, VW ft ) n , - (VV, VWfc)„« , VW h e W h (0). 
It is well known that VhV — > V strongly in Hl(Cl*). Given that e is uniformly bounded, 



n+l 



V n+1 )\\l < (e*( 
+ (e*( 



n+l' 



n+l 



b r h L+1 )W^ + \V {V h V n+1 - V n+1 )) 



>n+l\ 



m+1 



ft 

•n+l 



+<£*( < / , r i ) v ^ n+i ) v ( vr 
=i+n+ in. 



yn+l 



n+lj 

)>n« 



Let us estimate each one of the terms separately. Since the coefficients are bounded and the 
sequence VV^ +1 is uniformly bounded in L 2 (0), the strong convergence of VhV n+1 shows 
that I — > 0. For II we use the equation, namely 

--(Ql + \V£ +1 -V h V n+1 ) 



II 



<e*(^ +1 )V^ l+1 ,V(y, l " +1 



P h V n+1 )) 







since q h 
as 



1 converges strongly in L 2 (Q). Finally, notice that the last term can be rewritten 



111 = {{e*{(j} n h +1 )\ -e*{d? n+1 ) W n+1 ,V(V 



n+l 



yn+l 



Ln+1\ 



W" +i ,V(y 



n+l 



The uniform boundedness of VV^ +1 in L 2 (f2) implies that, for the first term, it suffices to 
show that (e*{(t>l +1 ) - e*((f> n+1 )) W n+1 -+ in L 2 (ft), which follows from the Lebesgue 
dominated convergence theorem. For the second term, use the weak convergence of VV^ +1 . 
This, together with the strong L 2 -convergence of q^ +1 implies that the limit solves (6.1). 
Equation (6.2): The strong L 2 -convergence of implies that — > ^Ug n+1 strongly 
in L 2 (f2). Using the compact embeddings <g L 4 (fi) and V <s L 4 (f2), we see that 

(« +1 ,Vr) -+ (<f u"+\ Vr) , Vr G ff 1 ^), 



ft 



n+l\ 



can be treated as in (6.8). These observations 



as h -> 0. The term iC(^)V(A g ; 
imply that the limit solves ( |6.2[ ). 
Equation (6.3): The strong L 2 -convergence of u^ +1 , the weak H 1 -convergence of <^ +1 and 
an argument similar to (6.8) imply that the limit solves (6.3). 
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Equation (6.4): The smoothness of W and the fact its growth is quadratic imply 
KW'(^) - W(<j> n ),P)\ < max|>V»|||^ - 0" 



0. 



A similar argument and the embedding d L 2 (T) can be used to show convergence 

of 9y s (<^). Since p is a bounded smooth function, 

(p'(<^K-<+\/2) -+ (p'(<TK-u"+\/2) . 
The strong L 2 -convergence of VV^" 1 " 1 implies that 

(£(<^ +1 ,<Ow,r +1 l 2 ^) -+ <f(0" +1 ,0")|VT/" +1 | 2 ,/2), 



where it is necessary to have fj, € L°°(f2). To conclude that (6.4) is satisfied by the limit, 



it is left to show that </>^ +1 converges strongly in L 2 (T). We know that <^ +1 converges 
weakly in L 2 (T). On the other hand -^f.~d4>h +1 converges strongly in £ 2 (T), uj^t converges 
strongly in L 2 (r) and ip((f>h) converges a.e. in T. 



Equations (6.5): Clearly, ( |6.5b ) is satisfied. To show that (6.5a) holds, notice that 



(p(^ +1 )< +1 - /3 (r +1 )u" +1 ,w) = 

(p(K +1 ) « +1 - u" +1 ) , w) + ( (p(^ +1 ) - PW^)) u n+1 , w) 0. 
Since we assume that tp is smooth and the slip coefficient (3 is smooth and depends only on 



the phase field, but not on the stress (as opposed to f 2.4 ) , we can get convergence of the 
terms [$(^)u^ , w T ] and [ujjj w x ^(^)] , respectively. The advection term can 

be treated using standard arguments and thus we will not give details here. The terms 

(^ +1 Vftw) , <tfV(A<# +1 + C^Vw) , 

can be treated using arguments similar to the ones given before. The term 



At 



can be easily shown to converge since all terms converge strongly. The convergence of the 
term 



At 



follows again from the compact embedding <<= L 2 (T). Finally, the convergence of the 

viscous stress term follows the lines of the proof of (|6. 



To conlcude, we notice that we do not need to reprove estimates similar to (4.16). These are 



uniformly valid, in h, for all terms in the sequence and, therefore, valid for the limit. Moreover, if 



one wanted to obtain an energy estimate by repeating the arguments used to obtain Proposition |4.14 
it would be necessary first to obtain uniform L°° bounds on the sequence D(j>h,At, since one of the 
steps in the proof requires setting p, h = 2d<j)^ +1 . □ 

Remark 6.9 (Limit At — > 0). We are not able to pass to the limit when At — > for several reasons. 
First, the estimates on the pressure depend on the time-step and getting around this would require 
finer estimates on the time derivative of the velocity, this is standard for Navier Stokes. In addition, 
the terms 

/ MP >..„ i 

\ At 



At 



-u , w 



-u 



,w 
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would require finer estimates on the time derivative of the phase field which we have not been able 
to show. It might be possible, however, to circumvent these two restrictions by defining the weak 
solution to the continuous problem with an unconstrained formulation for the momentum equation 
(i.e., solution and test functions in V) and modifying the Cahn-Hilliard equations to their "viscous 
version", in other words suitably adding a term of the form <p t . We will not pursue this direction. 

7. Conclusions and Perspectives 

Some possible directions for future work would be to extend the analysis by passing to the limit 
as At — > 0, or investigate the phenomenological pinning model more thoroughly. It would also be 
interesting to look at the use of open boundary conditions on <9*f2*, which is more physically correct 
for some electrowetting devices. As far as we know, this is an open area of research in the context of 
phase- field methods. Other extensions of the model could include the transport of surfactant at the 
liquid-gas interface, though this would make the model more complicated. We want to emphasize 
that our model gives physically reasonable results when modeling actual electrowetting systems, 
and so could be used within an optimization framework for improving device design. 

Concerning numerics, an important issue that has not been addressed is how to actually solve 
the discretized systems. Even in the fully uncoupled case, the pressence of the dynamic boundary 
condition in the Cahn-Hilliard system (Step 2 in the scheme of section [5| makes this problem 
extremely ill-conditioned and standard preconditioning techniques (for instance the one in [5]) 
inapplicable. 
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